home *** CD-ROM | disk | FTP | other *** search
/ PD Collection CD 1 / PD Collection CD 1.iso / programer2 / pari2 / pari / c / base < prev    next >
Text File  |  1991-12-06  |  36KB  |  1,430 lines
  1. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  2. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  3. /*                                                                 */
  4. /*                       BASE D'ENTIERS                            */
  5. /*                                                                 */
  6. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  7. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  8.  
  9. # include "genpari.h"
  10. GEN rquot(),ordmax(),rtran(),mtran(),matinv();
  11. void rowred();
  12.  
  13. GEN allbase(x,code,y)
  14.  
  15.      GEN x,*y;
  16.      long code;
  17.  
  18. /*******************************************************************
  19.   Entree:     x polynome unitaire a coefficients dans Z de deg n
  20.         definissant un corps de nombres K=Q(theta);
  21.               code 0, 1 ou (long)p selon que l'on veut base, smallbase
  22.             ou factoredbase.
  23.           y pointeur sur un GEN destine a recevoir
  24.         le discriminant du corps K.
  25.   Sortie:    retourne 1) un vecteur (horizontal) a n composantes, 
  26.             de polynomes a coeff dans Q (de deg 0,1...n-1)
  27.         constituant une base de l'anneau des entiers de K.
  28.             2) le discriminant de K (dans *y).
  29.         Rem: le denominateur commun des coef. est dans da.
  30. *******************************************************************/
  31.  
  32. {
  33.   GEN p,a,at,bt,b,da,db,q;
  34.   long av=avma,tetpil,n,h,j,je,k,r,s,t,pro,v;
  35.  
  36.   if(typ(x)!=10) err(allbaser1);
  37.   n=lgef(x)-3;if(n<=0) err(allbaser1);
  38.   v=varn(x);*y=discsr(x);
  39.   switch(code)
  40.     {
  41.     case 0: p=auxdecomp(absi(*y),1);h=lg(p[1])-1;break; /* base */
  42.     case 1: p=auxdecomp(absi(*y),0);h=lg(p[1])-1;break; /* smallbase */
  43.     default: p=(GEN)code;
  44.       if((typ(p)!=19)||(lg(p)!=3)) err(factoreder1); /* factoredbase */
  45.       h=lg(p[1])-1;
  46.       q=gun;for(je=1;je<=h;je++) q=gmul(q,gpui(coeff(p,je,1),coeff(p,je,2)));
  47.       if(gcmp(absi(q),absi(*y))) err(factoreder2);
  48.     }
  49.   a=idmat(n);da=gun;
  50.   for(je=1;je<=h;je++)
  51.     {
  52.       if(gcmpgs(coeff(p,je,2),1)>0)
  53.     {
  54.       b=ordmax(x,coeff(p,je,1),coeff(p,je,2),&db);
  55.       a=gmul(db,a);b=gmul(da,b);
  56.       da=mulii(db,da);db=da;
  57.       at=gtrans(a);bt=gtrans(b);
  58.       for(r=n-1;r>=0;r--)
  59.         for(s=r;s>=0;s--)
  60.           while(signe(coef1(bt,s,r)))
  61.         {
  62.           q=rquot(coef1(at,s,s),coef1(bt,s,r));
  63.           at[s+1]=(long)rtran(at[s+1],bt[r+1],q);
  64.           for(t=s-1;t>=0;t--)
  65.             {
  66.               q=rquot(coef1(at,t,s),coef1(at,t,t));
  67.               at[s+1]=(long)rtran(at[s+1],at[t+1],q);
  68.             }
  69.           pro=at[s+1];at[s+1]=bt[r+1];bt[r+1]=pro;
  70.         }
  71.       for(j=n-1;j>=0;j--)
  72.         {
  73.           for(k=0;k<j;k++)
  74.         {
  75.           while(signe(coef1(at,j,k)))
  76.             {
  77.               q=rquot(coef1(at,j,j),coef1(at,j,k));
  78.               at[j+1]=(long)rtran(at[j+1],at[k+1],q);
  79.               pro=at[j+1];at[j+1]=at[k+1];at[k+1]=pro;
  80.             }
  81.         }
  82.           if(signe(coef1(at,j,j))<0)
  83.         for(k=0;k<=j;k++) coef1(at,k,j)=lnegi(coef1(at,k,j));
  84.           for(k=j+1;k<n;k++)
  85.         {
  86.           q=rquot(coef1(at,j,k),coef1(at,j,j));
  87.           at[k+1]=(long)rtran(at[k+1],at[j+1],q);
  88.         }
  89.         }
  90.       for(j=1;j<n;j++)
  91.         if(!cmpii(coef1(at,j,j),coef1(at,j-1,j-1)))
  92.           {
  93.         coef1(at,0,j)=zero;
  94.         for(k=1;k<=j;k++)
  95.           coef1(at,k,j)=coef1(at,k-1,j-1);
  96.           }
  97.       a=gtrans(at);
  98.     }
  99.     }
  100.   for(j=1;j<=n;j++)
  101.     {
  102.       *y=divii(mulii(coeff(a,j,j),*y),da);
  103.       *y=divii(mulii(coeff(a,j,j),*y),da);
  104.     }
  105.   tetpil=avma;*y=gcopy(*y);at=cgetg(n+1,17);
  106.   for(j=1;j<=n;j++)
  107.     {
  108.       q=cgetg(j+2,10);q[1]=0x1000002+j+(v<<16);at[j]=(long)q;
  109.       for(k=2;k<=j+1;k++) q[k]=ldiv(coeff(a,j,k-1),da);
  110.     }
  111.   pro=lpile(av,tetpil,0)>>2;at+=pro;(*y)+=pro;
  112.   return at;
  113. }
  114.  
  115. GEN base(x,y)
  116.      GEN x,*y;
  117. {
  118.   return allbase(x,0,y);
  119. }
  120.  
  121. GEN smallbase(x,y)
  122.      GEN x,*y;
  123. {
  124.   return allbase(x,1,y);
  125. }
  126.  
  127. GEN factoredbase(x,p,y)
  128.      GEN x,p,*y;
  129. {
  130.   return allbase(x,(long)p,y);
  131. }
  132.  
  133. GEN discf(x)
  134.      GEN x;
  135. {
  136.   GEN y;
  137.   long av,tetpil;
  138.  
  139.   av=avma;allbase(x,0,&y);tetpil=avma;
  140.   return gerepile(av,tetpil,gcopy(y));
  141. }
  142.  
  143. GEN smalldiscf(x)
  144.      GEN x;
  145. {
  146.   GEN y;
  147.   long av,tetpil;
  148.  
  149.   av=avma;allbase(x,1,&y);tetpil=avma;
  150.   return gerepile(av,tetpil,gcopy(y));
  151. }
  152.  
  153. GEN factoreddiscf(x,p)
  154.      GEN x,p;
  155. {
  156.   GEN y;
  157.   long av,tetpil;
  158.  
  159.   av=avma;allbase(x,(long)p,&y);tetpil=avma;
  160.   return gerepile(av,tetpil,gcopy(y));
  161. }
  162.  
  163. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  164. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  165. /*                                                                 */
  166. /*   Quotient et Reste normalises   ( -1/2 < r = x-q*y <= 1/2 )    */
  167. /*                                                                 */
  168. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  169. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  170.  
  171. GEN rquot(x,y)
  172.  
  173.      GEN x,y;
  174.  
  175. {
  176.   GEN u,v,w,p;
  177.   long av,av1;
  178.  
  179.   av=avma;
  180.   u=absi(y);v=shifti(x,1);w=shifti(y,1);
  181.   if ( cmpii(u,v)>0) p=subii(v,u);
  182.   else p=addsi(-1,addii(u,v));
  183.   av1=avma;
  184.   return(gerepile(av,av1,divii(p,w)));
  185. }
  186.  
  187. GEN rrmdr(x,y)
  188.  
  189.      GEN x,y;
  190.  
  191. {
  192.   GEN p;
  193.   long av,av1;
  194.  
  195.   av=avma;
  196.   p=mulii(rquot(x,y),y);
  197.   av1=avma;
  198.   return(gerepile(av,av1,subii(x,p)));
  199. }
  200.  
  201. GEN rinv(x,y)
  202.  
  203.      GEN x,y;
  204.  
  205. {
  206.   GEN a,c,q,r,t;
  207.   long av,av1;
  208.  
  209.   av=avma;
  210.   a=gun;c=gzero;
  211.   while( signe(y))
  212.     {
  213.       q=rquot(x,y);
  214.       r=subii(a,mulii(q,c));a=c;c=r;
  215.       t=subii(x,mulii(q,y));x=y;y=t;
  216.     }
  217.   av1=avma;
  218.   if (signe(x)<0) a=negi(a);
  219.   if (signe(c)) { av1=avma; a=rrmdr(a,c); }
  220.   return( gerepile(av,av1,a));
  221. }
  222.  
  223. GEN rgcd(x,y)
  224.  
  225.      GEN x,y;
  226.  
  227. {
  228.   GEN z;
  229.   long av,av1;
  230.  
  231.   av=avma;
  232.   while(signe(y))
  233.     {
  234.       z=rrmdr(x,y);x=y;y=z;
  235.     }
  236.   av1=avma;
  237.   return(gerepile(av,av1,absi(x)));
  238. }
  239.     
  240.  
  241. GEN rlcm(x,y)
  242.  
  243.      GEN x,y;
  244.  
  245. {
  246.   GEN d,z;
  247.   long av,av1;
  248.  
  249.   av=avma;
  250.   z=mulii(x,y);d=rgcd(x,y);
  251.   av1=avma;
  252.   return(gerepile(av,av1,divii(z,d)));
  253. }
  254.  
  255. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  256. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  257. /*                                                                 */
  258. /*           Matrice compagnon du polynome unitaire x              */
  259. /*                                                                 */
  260. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  261. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  262.  
  263. GEN companion(x)
  264.  
  265.      GEN     x;
  266.  
  267. {
  268.   long    i,j,l;
  269.   GEN     y;
  270.   
  271.   l=lgef(x)-2;y=cgetg(l,19);
  272.   for(i=1;i<l;i++) y[i]=lgetg(l,18);
  273.   for(i=0;i<l-2;i++)
  274.     for(j=0;j<l-1;j++) coef1(y,i,j)=((i+1)==j) ? un : zero;
  275.   for(j=0;j<l-1;j++) coef1(y,l-2,j)=lneg(x[j+2]);
  276.   return(y);
  277. }
  278.  
  279.  
  280.  
  281. GEN ordmax(f,p,e,ptdelta)
  282.  
  283. GEN f,p,e;
  284. GEN *ptdelta;
  285.  
  286. {
  287.   GEN m,hh,pp,dd,ppdd,index,q,r,s,b,c,t,jp,v,delta;
  288.   GEN cf[20],w[20],a;
  289.   long h,i,j,k,sp,epsilon,n=lgef(f)-3;
  290.  
  291.   a=cgetg(n*n+1,19);
  292.   for(j=1;j<=n*n;j++)
  293.   {
  294.     a[j]=lgetg(n+1,18);
  295.     for(k=1;k<=n;k++) coeff(a,k,j)=zero;
  296.   }
  297.   v=cgetg(n+1,18);
  298.   cf[0]=idmat(n);
  299.   cf[1]=companion(f);
  300.   for(j=2;j<n;j++) cf[j]=gmul(cf[1],cf[j-1]);
  301.   delta=gun; epsilon=itos(e);
  302.   m=idmat(n);
  303.  
  304.   do
  305.     {
  306.       pp=mulii(p,p);
  307.       dd=mulii(delta,delta);
  308.       ppdd=mulii(dd,pp);
  309.       b=matinv(m,delta);
  310.       for(i=0;i<n;i++)
  311.     {
  312.       t=gscalsmat(0,n); /* t <--- matrice nulle d'ordre n */
  313.       for(h=0;h<n;h++)
  314.         for(j=0;j<n;j++)
  315.           for(k=0;k<n;k++)
  316.         coef1(t,j,k)=(long)rrmdr(addii(coef1(t,j,k),mulii(coef1(m,i,h),coef1(cf[h],j,k))),ppdd);
  317.       c=gmul(t,b);
  318.       w[i]=gmul(m,c);
  319.       for(j=0;j<n;j++)
  320.         for(k=0;k<n;k++)
  321.           coef1(w[i],j,k)=(long)rrmdr(divii(coef1(w[i],j,k),dd),pp);
  322.     }
  323.       if(cmpis(p,n)>0)
  324.     {
  325.       for(i=0;i<n;i++)
  326.         for(j=0;j<n;j++)
  327.           {
  328.         coeff(t,i+1,j+1)=zero;
  329.         for(k=0;k<n;k++)
  330.           for(h=0;h<n;h++)
  331.             {
  332.               r=modii(coef1(w[i],k,h),p);
  333.               s=modii(coef1(w[j],h,k),p);
  334.               coef1(t,i,j)=lmodii(addii(coef1(t,i,j),mulii(r,s)),p);
  335.             }
  336.           }
  337.     }
  338.       else
  339.     {
  340.       for(j=0;j<n;j++)
  341.         {
  342.           for(i=0;i<n;i++)
  343.         coef1(b,i,j)=(i==0)? un:zero;
  344. /* ici la boucle en k calcule la puissance p mod p de w[j] */
  345.           sp=itos(p);
  346.           for(k=0;k<sp;k++)
  347.         {
  348.           for(i=0;i<n;i++)
  349.             {
  350.               v[i+1]=zero;
  351.               for(h=0;h<n;h++)
  352.             v[i+1]=lmodii(addii(v[i+1],mulii(coef1(b,h,j),coef1(w[j],h,i))),p);
  353.             }
  354.           for(i=0;i<n;i++) coef1(b,i,j)=v[i+1];
  355.         }
  356.         }
  357.       q=p;t=b;
  358.       while(cmpis(q,n)<0)
  359.         {
  360.           q=mulii(q,p);
  361.           t=gmul(b,t);
  362.         }
  363.     }
  364.       for(i=0;i<n;i++)
  365.     for(j=0;j<n;j++)
  366.       {
  367.         coef1(a,j,i)=(i==j)? (long)p:zero;
  368.         coef1(a,j,n+i)=lmodii(coef1(t,i,j),p);
  369.       }
  370.       rowred(a,2*n-1,pp);
  371.       for(i=0;i<n;i++)
  372.     for(j=0;j<n;j++)
  373.       coef1(b,j,i)=coef1(a,j,i);
  374.       jp=matinv(b,p);
  375.       for(k=0;k<n;k++)
  376.     {
  377.       t=gmul(jp,w[k]);
  378.       t=gmul(t,b);
  379.       for(i=0;i<n;i++)
  380.         for(j=0;j<n;j++)
  381.           coef1(t,i,j)=ldivii(coef1(t,i,j),p);
  382.       h=0;
  383.       for(i=0;i<n;i++)
  384.         for(j=0;j<n;j++)
  385.           {
  386.         coef1(a,k,h)=coef1(t,i,j);
  387.         h++;
  388.           }
  389.     }
  390.       rowred(a,n*n-1,pp);
  391.       index=gun;
  392.       for(i=0;i<n;i++)
  393.     index=mulii(index,coef1(a,i,i));
  394.       if (cmpsi(1,index))
  395.     {
  396.       delta=mulii(index,delta);
  397.       for(i=0;i<n;i++)
  398.         for(j=0;j<n;j++)
  399.           coef1(c,i,j)=coef1(a,i,j);
  400.       b=matinv(c,index);
  401.       m=gmul(b,m);
  402.       hh=delta;
  403.       for(i=0;i<n;i++)
  404.         for(j=0;j<n;j++)
  405.           hh=rgcd(coef1(m,i,j),hh);
  406.       if(cmpis(hh,1)>1)
  407.         {
  408.           delta=divii(delta,hh);
  409.           for(i=0;i<n;i++)
  410.         for(j=0;j<n;j++)
  411.           coef1(m,i,j)=ldivii(coef1(m,i,j),hh);
  412.         }
  413.       q=index;
  414.       while(!signe(modii(q,p)))
  415.         {
  416.           q=divii(q,p);
  417.           epsilon=epsilon-2;
  418.         }
  419.     }
  420.     }
  421.   while(!gcmp1(index) && (epsilon>=2));
  422.   *ptdelta=delta;
  423.   return(m);
  424. }
  425.  
  426. void rowred(a,rlim,rmod)
  427.      GEN a,rmod;
  428.      long rlim;
  429.  
  430. {
  431.   long j,k,n,pro;
  432.   GEN q;
  433.  
  434.   n=lg(a[1])-1;
  435.   for(j=0;j<n;j++)
  436.     {
  437.       for(k=j+1;k<=rlim;k++)
  438.     while (signe(coef1(a,j,k)))
  439.       {
  440.         q=rquot(coef1(a,j,j),coef1(a,j,k));
  441.         a[j+1]=(long)mtran(a[j+1],a[k+1],q,rmod);
  442.         pro=a[j+1];a[j+1]=a[k+1];a[k+1]=pro;
  443.       }
  444.       if (signe(coef1(a,j,j))<0)
  445.     for(k=j;k<n;k++) coef1(a,k,j)=lnegi(coef1(a,k,j));
  446.       for(k=0;k<j;k++)
  447.     {
  448.       q=rquot(coef1(a,j,k),coef1(a,j,j));
  449.       a[k+1]=(long)mtran(a[k+1],a[j+1],q,rmod);
  450.     }
  451.     }
  452. }
  453.  
  454. GEN rtran(v,w,q)
  455.      GEN v,w,q ;
  456.  
  457. {
  458.   long av,tetpil;
  459.   GEN p1;
  460.  
  461.   if (signe(q))
  462.     {
  463.       av=avma;p1=gmul(q,w);tetpil=avma;
  464.       return gerepile(av,tetpil,gsub(v,p1));
  465.     }
  466.   else return v;
  467. }
  468.  
  469. GEN mtran(v,w,q,m)
  470.      GEN v,w,q,m;
  471.  
  472. {
  473.   long k;
  474.   
  475.   if (signe(q))
  476.     {
  477.       for(k=0;k<lg(v)-1;k++)
  478.     {
  479.       v[k+1]=(long)rrmdr(subii(v[k+1],modii(mulii(q,w[k+1]),m)),m);
  480.     }
  481.     }
  482.   return v;
  483. }
  484.  
  485.  
  486. GEN matinv(x,d)
  487.      GEN x,d;
  488. /*=======================================================================
  489.     Calcule d/x  ou  d est entier et x matrice triangulaire inferieure
  490.   entiere dont les coeff diagonaux divisent d ( resultat entier).
  491. ========================================================================*/
  492. {
  493.   long n,h,i,j,k,av,av1;
  494.   GEN y;
  495.  
  496.   av=avma;
  497.   y=idmat(n=lg(x)-1);
  498.   for(i=1;i<=n;i++)
  499.     coeff(y,i,i)=ldivii(d,coeff(x,i,i));
  500.   for(i=2;i<=n;i++)
  501.     for(j=i-1;j;j--)
  502.       {
  503.     for(h=zero,k=j+1;k<=i;k++)
  504.       h=ladd(h,mulii(coeff(y,i,k),coeff(x,k,j)));
  505.     coeff(y,i,j)=ldivii(negi(h),coeff(x,j,j));
  506.       }
  507.   av1=avma;
  508.   return gerepile(av,av1,gcopy(y));
  509. }
  510.  
  511. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  512. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  513. /*~                                    ~*/
  514. /*~            HERMITE NORMAL FORM REDUCTION            ~*/
  515. /*~                                    ~*/
  516. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  517. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  518.  
  519. GEN hnfold(x)
  520.      GEN x;
  521. {
  522.   long li,co,av,tetpil,i,j,j1,def,ind,lim;
  523.   GEN p1,p2,y,z,m1;
  524.  
  525.   if(typ(x)!=19) err(hnfer1);
  526.   lim=(avma+bot)>>1;
  527.   av=avma;co=lg(x);li=lg(x[1]);y=gcopy(x);
  528.   def=co;
  529.   if(li>co) err(hnfer2);
  530.   for(i=li-1;i>=1;i--)
  531.     {
  532.       def--;j=1;while((j<def)&&(!signe(coeff(y,i,j)))) j++;
  533.       while(j<def)
  534.     {
  535.       m1=absi(coeff(y,i,j));ind=j;
  536.       for(j1=j+1;j1<def;j1++)
  537.         {
  538.           p1=(GEN)coeff(y,i,j1);
  539.           if(signe(p1)&&(cmpii(absi(p1),m1)<0))
  540.         {
  541.           m1=p1;ind=j1;
  542.         }
  543.         }
  544.       p1=(GEN)coeff(y,i,def);
  545.       if(!(signe(p1)&&(cmpii(absi(p1),m1)<0)))
  546.         {
  547.           p1=(GEN)y[def];y[def]=y[ind];y[ind]=(long)p1;
  548.         }
  549.       p1=(GEN)coeff(y,i,def);if(signe(p1)<0) y[def]=lneg(y[def]);
  550.       p1=(GEN)coeff(y,i,def);
  551.       for(j=1;j<def;j++)
  552.         {
  553.           p2=gdivent(coeff(y,i,j),p1);
  554.           y[j]=lsub(y[j],gmul(p2,y[def]));
  555.         }
  556.       j=1;while((j<def)&&(!signe(coeff(y,i,j)))) j++;
  557.     }
  558.       p1=(GEN)coeff(y,i,def);
  559.       if(signe(p1)<0) {y[def]=lneg(y[def]);p1=(GEN)coeff(y,i,def);}
  560.       if(signe(p1))
  561.     {
  562.       for(j=def+1;j<co;j++)
  563.         {
  564.           p2=gdivent(coeff(y,i,j),p1);
  565.           y[j]=lsub(y[j],gmul(p2,y[def]));
  566.         }
  567.     }
  568.       if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
  569.     }
  570.   tetpil=avma;z=cgetg(li,19);
  571.   for(j=1;j<li;j++)
  572.     {
  573.       z[j]=lcopy(y[j+co-li]);
  574.     }
  575.   for(j=1;j<=co-li;j++) if(!gcmp0(y[j])) err(hnfer3);
  576.   return gerepile(av,tetpil,z);
  577. }
  578.  
  579. GEN hnfold2(x)
  580.      GEN x;
  581. {
  582.   long li,co,av,tetpil,i,j,def,lim;
  583.   GEN p1,p2,y,z,u,v,d;
  584.  
  585.   if(typ(x)!=19) err(hnfer1);
  586.   lim=(avma+bot)>>1;
  587.   av=avma;co=lg(x);li=lg(x[1]);y=gcopy(x);
  588.   def=co;
  589.   if(li>co) err(hnfer2);
  590.   for(i=li-1;i>=1;i--)
  591.     {
  592.       def--;j=def;while(j&&(!signe(coeff(y,i,j)))) j--;
  593.       while(j>1)
  594.     {
  595.       d=bezout(coeff(y,i,j),coeff(y,i,j-1),&u,&v);
  596.       p1=gadd(gmul(u,y[j]),gmul(v,y[j-1]));
  597.       y[j]=lsub(gmul(divii(coeff(y,i,j),d),y[j-1]),gmul(divii(coeff(y,i,j-1),d),y[j]));
  598.       y[j-1]=(long)p1;
  599.       j--;while(j&&(!signe(coeff(y,i,j)))) j--;
  600.     }
  601.       if(j==1)
  602.     {
  603.       d=bezout(coeff(y,i,1),coeff(y,i,def),&u,&v);
  604.       p1=gadd(gmul(u,y[1]),gmul(v,y[def]));
  605.       y[1]=lsub(gmul(divii(coeff(y,i,1),d),y[def]),gmul(divii(coeff(y,i,def),d),y[1]));
  606.       y[def]=(long)p1;
  607.     }
  608.       p1=(GEN)coeff(y,i,def);
  609.       if(signe(p1)<0) {y[def]=lneg(y[def]);p1=(GEN)coeff(y,i,def);}
  610.       if(signe(p1))
  611.     {
  612.       for(j=def+1;j<co;j++)
  613.         {
  614.           p2=gdivent(coeff(y,i,j),p1);
  615.           y[j]=lsub(y[j],gmul(p2,y[def]));
  616.         }
  617.     }
  618.       if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
  619.     }
  620.   tetpil=avma;z=cgetg(li,19);
  621.   for(j=1;j<li;j++)
  622.     {
  623.       z[j]=lcopy(y[j+co-li]);
  624.     }
  625.   for(j=1;j<=co-li;j++) if(!gcmp0(y[j])) err(hnfer3);
  626.   return gerepile(av,tetpil,z);
  627. }
  628.  
  629. GEN hnf(x)
  630.      GEN x;
  631. {
  632.   long li,co,av,tetpil,i,j,k,def,ldef,lim;
  633.   GEN p1,p2,y,z,u,v,d;
  634.  
  635.   if(typ(x)!=19) err(hnfer1);
  636.   lim=(avma+bot)>>1;
  637.   av=avma;co=lg(x);li=lg(x[1]);y=gcopy(x);
  638.   def=co;ldef=(li>co)?li-co+1:1;
  639.   for(i=li-1;i>=ldef;i--)
  640.     {
  641.       def--;j=def-1;while(j&&(!signe(coeff(y,i,j)))) j--;
  642.       while(j>1)
  643.     {
  644.       d=bezout(coeff(y,i,j),coeff(y,i,j-1),&u,&v);
  645.       p1=gadd(gmul(u,y[j]),gmul(v,y[j-1]));
  646.       y[j]=lsub(gmul(divii(coeff(y,i,j),d),y[j-1]),gmul(divii(coeff(y,i,j-1),d),y[j]));
  647.       y[j-1]=(long)p1;
  648.       j--;while(j&&(!signe(coeff(y,i,j)))) j--;
  649.     }
  650.       if(j==1)
  651.     {
  652.       d=bezout(coeff(y,i,1),coeff(y,i,def),&u,&v);
  653.       p1=gadd(gmul(u,y[1]),gmul(v,y[def]));
  654.       y[1]=lsub(gmul(divii(coeff(y,i,1),d),y[def]),gmul(divii(coeff(y,i,def),d),y[1]));
  655.       y[def]=(long)p1;
  656.     }
  657.       p1=(GEN)coeff(y,i,def);
  658.       if(signe(p1)<0) {y[def]=lneg(y[def]);p1=(GEN)coeff(y,i,def);}
  659.       if(signe(p1))
  660.     {
  661.       for(j=def+1;j<co;j++)
  662.         {
  663.           p2=gdivent(coeff(y,i,j),p1);
  664.           y[j]=lsub(y[j],gmul(p2,y[def]));
  665.         }
  666.     }
  667.       else def++;
  668.       if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
  669.     }
  670.   for(i=0,j=1;j<co;j++) if(!gcmp0(y[j])) i++;
  671.   tetpil=avma;z=cgetg(i+1,19);
  672.   for(k=0,j=1;j<co;j++) if(!gcmp0(y[j])) z[++k]=lcopy(y[j]);
  673.   return gerepile(av,tetpil,z);
  674. }
  675.  
  676.  
  677. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  678. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  679. /*~                                    ~*/
  680. /*~            SMITH NORMAL FORM REDUCTION            ~*/
  681. /*~                                    ~*/
  682. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  683. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  684.  
  685. GEN smith2(x)
  686.      GEN x;
  687. {
  688.   long li,av,tetpil,i,j,k,l,lim,c,fl,n;
  689.   GEN p1,p2,p3,p4,y,z,b,u,v,d;
  690.  
  691.   if(typ(x)!=19) err(hnfer1);
  692.   lim=(avma+bot)>>1;
  693.   av=avma;n=lg(x)-1;li=lg(x[1])-1;y=gcopy(x);
  694.   if(li!=n) err(hnfer2);
  695.   for(i=n;i>=2;i--)
  696.     {
  697.       do
  698.     {
  699.       for(j=i-1;j>=1;j--)
  700.         {
  701.           p1=(GEN)coeff(y,i,j);
  702.           if(signe(p1))
  703.         {
  704.           p2=(GEN)coeff(y,i,i);
  705.           if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
  706.           else if(!signe(addii(p1,p2))) 
  707.             {
  708.               d=absi(p1);u=(signe(p2)>0)?gun:gneg(un);
  709.               v=gzero;p3=u;p4=gneg(u);
  710.             }
  711.           else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
  712.           for(k=1;k<=i;k++)
  713.             {
  714.               b=addii(mulii(u,coeff(y,k,i)),mulii(v,coeff(y,k,j)));
  715.               coeff(y,k,j)=lsubii(mulii(p3,coeff(y,k,j)),mulii(p4,coeff(y,k,i)));
  716.               coeff(y,k,i)=(long)b;
  717.             }
  718.         }
  719.         }
  720.       c=0;
  721.       for(j=i-1;j>=1;j--)
  722.         {
  723.           p1=(GEN)coeff(y,j,i);
  724.           if(signe(p1))
  725.         {
  726.           p2=(GEN)coeff(y,i,i);
  727.           if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
  728.           else if(!signe(addii(p1,p2))) 
  729.             {
  730.               d=absi(p1);u=(signe(p2)>0)?gun:gneg(un);
  731.               v=gzero;p3=u;p4=gneg(u);
  732.             }
  733.           else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
  734.           for(k=1;k<=i;k++)
  735.             {
  736.               b=addii(mulii(u,coeff(y,i,k)),mulii(v,coeff(y,j,k)));
  737.               coeff(y,j,k)=lsubii(mulii(p3,coeff(y,j,k)),mulii(p4,coeff(y,i,k)));
  738.               coeff(y,i,k)=(long)b;
  739.             }
  740.           c++;
  741.         }
  742.         }
  743.       if(!c)
  744.         {
  745.           b=(GEN)coeff(y,i,i);fl=1;
  746.           if(signe(b))
  747.         {
  748.           for(k=1;(k<i)&&fl;k++)
  749.             for(l=1;(l<i)&&fl;l++)
  750.               fl= !signe(modii(coeff(y,k,l),b));
  751.         }
  752.           if(!fl) {l--;y[i]=ladd(y[i],y[l]);}
  753.         }
  754.       if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
  755.     }
  756.       while(c||(!fl));
  757.     }
  758.   tetpil=avma;z=cgetg(n+1,17);
  759.   for(k=1;k<=n;k++) z[k]=labs(coeff(y,k,k));
  760.   return gerepile(av,tetpil,z);
  761. }
  762.  
  763. GEN smith(x)
  764.      GEN x;
  765. {
  766.   long li,av,tetpil,i,j,k,l,lim,c,fl,n;
  767.   GEN p1,p2,p3,p4,y,z,b,u,v,d;
  768.  
  769.   if(typ(x)!=19) err(hnfer1);
  770.   lim=(avma+bot)>>1;
  771.   av=avma;n=lg(x)-1;if(!n) return cgetg(1,17);
  772.   li=lg(x[1])-1;y=gcopy(x);
  773.   if(li!=n) err(hnfer2);
  774.   for(i=n;i>=2;i--)
  775.     {
  776.       do
  777.     {
  778.       c=0;
  779.       for(j=i-1;j>=1;j--)
  780.         {
  781.           p1=(GEN)coeff(y,i,j);
  782.           if(signe(p1))
  783.         {
  784.           p2=(GEN)coeff(y,i,i);
  785.           if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
  786.           else if(!signe(addii(p1,p2))) 
  787.             {
  788.               d=absi(p1);u=(signe(p2)>0)?gun:gneg(un);
  789.               v=gzero;p3=u;p4=gneg(u);
  790.             }
  791.           else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
  792.           for(k=1;k<=i;k++)
  793.             {
  794.               b=addii(mulii(u,coeff(y,k,i)),mulii(v,coeff(y,k,j)));
  795.               coeff(y,k,j)=lsubii(mulii(p3,coeff(y,k,j)),mulii(p4,coeff(y,k,i)));
  796.               coeff(y,k,i)=(long)b;
  797.             }
  798.           c++;
  799.         }
  800.         }
  801.       for(j=i-1;j>=1;j--)
  802.         {
  803.           p1=(GEN)coeff(y,j,i);
  804.           if(signe(p1))
  805.         {
  806.           p2=(GEN)coeff(y,i,i);
  807.           if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
  808.           else if(!signe(addii(p1,p2))) 
  809.             {
  810.               d=absi(p1);u=(signe(p2)>0)?gun:gneg(un);
  811.               v=gzero;p3=u;p4=gneg(u);
  812.             }
  813.           else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
  814.           for(k=1;k<=i;k++)
  815.             {
  816.               b=addii(mulii(u,coeff(y,i,k)),mulii(v,coeff(y,j,k)));
  817.               coeff(y,j,k)=lsubii(mulii(p3,coeff(y,j,k)),mulii(p4,coeff(y,i,k)));
  818.               coeff(y,i,k)=(long)b;
  819.             }
  820.         }
  821.         }
  822.       if(!c)
  823.         {
  824.           b=(GEN)coeff(y,i,i);fl=1;
  825.           if(signe(b))
  826.         {
  827.           for(k=1;(k<i)&&fl;k++)
  828.             for(l=1;(l<i)&&fl;l++)
  829.               fl= !signe(modii(coeff(y,k,l),b));
  830.         }
  831.           if(!fl)
  832.         {
  833.           k--;
  834.           for(l=1;l<=i;l++)
  835.             coeff(y,i,l)=laddii(coeff(y,i,l),coeff(y,k,l));
  836.         }
  837.         }
  838.       if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
  839.     }
  840.       while(c||(!fl));
  841.     }
  842.   tetpil=avma;z=cgetg(n+1,17);
  843.   for(j=0,k=1;k<=n;k++) if(!signe(coeff(y,k,k))) z[++j]=zero;
  844.   for(k=1;k<=n;k++) if(signe(coeff(y,k,k))) z[++j]=labs(coeff(y,k,k));
  845.   return gerepile(av,tetpil,z);
  846. }
  847.  
  848.  
  849. GEN oldsmith(x)
  850.      GEN x;
  851. {
  852.   long li,av,tetpil,i,j,k,l,lim,c,fl,n;
  853.   GEN p1,p2,p3,p4,y,z,b,u,v,d;
  854.  
  855.   if(typ(x)!=19) err(hnfer1);
  856.   lim=(avma+bot)>>1;
  857.   av=avma;n=lg(x)-1;li=lg(x[1])-1;y=gcopy(x);
  858.   if(li!=n) err(hnfer2);
  859.   for(i=n;i>=2;i--)
  860.     {
  861.       do
  862.     {
  863.       c=0;
  864.       for(j=i-1;j>=1;j--)
  865.         {
  866.           p1=(GEN)coeff(y,i,j);
  867.           if(signe(p1))
  868.         {
  869.           p2=(GEN)coeff(y,i,i);
  870.           if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
  871.           else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
  872.           for(k=1;k<=i;k++)
  873.             {
  874.               b=addii(mulii(u,coeff(y,k,i)),mulii(v,coeff(y,k,j)));
  875.               coeff(y,k,j)=lsubii(mulii(p3,coeff(y,k,j)),mulii(p4,coeff(y,k,i)));
  876.               coeff(y,k,i)=(long)b;
  877.             }
  878.           c++;
  879.         }
  880.         }
  881.       for(j=i-1;j>=1;j--)
  882.         {
  883.           p1=(GEN)coeff(y,j,i);
  884.           if(signe(p1))
  885.         {
  886.           p2=(GEN)coeff(y,i,i);
  887.           if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
  888.           else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
  889.           for(k=1;k<=i;k++)
  890.             {
  891.               b=addii(mulii(u,coeff(y,i,k)),mulii(v,coeff(y,j,k)));
  892.               coeff(y,j,k)=lsubii(mulii(p3,coeff(y,j,k)),mulii(p4,coeff(y,i,k)));
  893.               coeff(y,i,k)=(long)b;
  894.             }
  895.         }
  896.         }
  897.       if(!c)
  898.         {
  899.           b=(GEN)coeff(y,i,i);fl=1;
  900.           if(signe(b))
  901.         {
  902.           for(k=1;(k<i)&&fl;k++)
  903.             for(l=1;(l<i)&&fl;l++)
  904.               fl= !signe(modii(coeff(y,k,l),b));
  905.         }
  906.           if(!fl)
  907.         {
  908.           k--;
  909.           for(l=1;l<=i;l++)
  910.             coeff(y,i,l)=laddii(coeff(y,i,l),coeff(y,k,l));
  911.         }
  912.         }
  913.       if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
  914.     }
  915.       while(c||(!fl));
  916.     }
  917.   tetpil=avma;z=cgetg(n+1,17);
  918.   for(j=0,k=1;k<=n;k++) if(!signe(coeff(y,k,k))) z[++j]=zero;
  919.   for(k=1;k<=n;k++) if(signe(coeff(y,k,k))) z[++j]=labs(coeff(y,k,k));
  920.   return gerepile(av,tetpil,z);
  921. }
  922.  
  923.  
  924. GEN transroot(x,i,j)
  925.      GEN x;
  926.      int i,j;
  927. {
  928.   long n=lg(x),k;
  929.   GEN y;
  930.  
  931.   y=cgetg(n,18);
  932.   for(k=1;k<n;k++) y[k]=((k==i)||(k==j))?x[i+j-k]:x[k];
  933.   return y;
  934. }
  935.  
  936. GEN tchirnhausen(x)
  937.      GEN x;
  938. {
  939.   long av=avma,tetpil,v,n,a,b,c;
  940.   GEN u;
  941.  
  942.   if(typ(x)!=10) err(galer1);
  943.   n=lgef(x)-3;if(n<=0) err(galer1);v=varn(x);
  944.   if(v) {u=gcopy(x);setvarn(u,0);x=u;}
  945.   do
  946.     {
  947.       a=rand()&15;if(!a) a=1;b=rand()&31;if(b>=16) b-=32;
  948.       c=rand()&31;if(c>=16) c-=32;
  949.       u=caract(gmodulcp(gaddsg(c,gmul(polx[0],gaddsg(b,gmulsg(a,polx[0])))),x),v);
  950.     }
  951.   while(lgef(polgcd(u,deriv(u,v)))>=4);
  952.   tetpil=avma;return gerepile(av,tetpil,gcopy(u));
  953. }
  954.  
  955. GEN galois(x,prec)
  956.      GEN x;
  957.      long prec;
  958. {
  959.  
  960.   long av=avma,av1,i,j,k,n,f,l,l2,e,e1;
  961.   GEN x1,p1,p2,p3,p4,p5,p6,y,z;
  962.   static int ind5[20]={2,5,3,4,1,3,4,5,1,5,2,4,1,2,3,5,1,4,2,3};
  963.   static int ind6[60]={3,5,4,6,2,6,4,5,2,3,5,6,2,4,3,6,2,5,3,4,1,4,5,6,1,5,3,6,1,6,3,4,1,3,4,5,1,6,2,5,1,2,4,6,1,5,2,4,1,3,2,6,1,2,3,5,1,4,2,3};
  964.  
  965.   if(typ(x)!=10) err(galer1);
  966.   n=lgef(x)-3;if(n<=0) err(galer1);
  967.   if(n>7) err(impl,"galois of degree higher than 7");
  968.   x=gdiv(x,content(x));
  969.   for(i=2;i<=n+2;i++) if(typ(x[i])!=1) err(galer2);
  970.   p1=(GEN)x[n+2];
  971.   if(!gcmp1(p1))
  972.     {
  973.       x1=cgetg(n+3,10);x1[1]=x[1];x1[n+2]=un;p2=gun;
  974.       for(i=n+1;i>=2;i--) {x1[i]=lmul(x[i],p2);if(i>2) p2=gmul(p1,p2);}
  975.       x=x1;
  976.     }
  977.   p1=factor(x);
  978.   if((lg(p1[1])>2)||(!gcmp1(coeff(p1,1,2)))) err(impl,"galois of reducible polynomial");
  979.   x1=gcopy(x);av1=avma;
  980.   for(;;)
  981.     {
  982.       switch(n)
  983.     {
  984.     case 1: avma=av;y=cgetg(4,17);y[1]=y[3]=un;y[2]=lneg(gun);return y;
  985.     case 2: avma=av;y=cgetg(4,17);y[1]=deux;y[3]=un;y[2]=lneg(gun);return y;
  986.     case 3: f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
  987.       if(f) {y[2]=un;y[1]=lstoi(3);return y;}
  988.       else {y[2]=lneg(gun);y[1]=lstoi(6);return y;}
  989.     case 4: 
  990.       do
  991.         {
  992.           p1=roots(x,prec);p2=p1;
  993.           p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
  994.           p4=gsub(polx[0],p3);p2=transroot(p1,1,2);
  995.           p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
  996.           p4=gmul(p4,gsub(polx[0],p3));p2=transroot(p1,1,3);
  997.           p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
  998.           p4=gmul(p4,gsub(polx[0],p3));p2=transroot(p1,1,4);
  999.           p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
  1000.           p4=gmul(p4,gsub(polx[0],p3));p2=transroot(p1,2,3);
  1001.           p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
  1002.           p4=gmul(p4,gsub(polx[0],p3));p2=transroot(p1,3,4);
  1003.           p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
  1004.           p4=gmul(p4,gsub(polx[0],p3));p5=grndtoi(greal(p4),&e);
  1005.           e=max(e,gexpo(gimag(p4)));
  1006.           if(e> -10) prec=(prec<<2)-2;
  1007.         }
  1008.       while(e> -10);
  1009.       p6=ggcd(p5,deriv(p5,0));
  1010.       f=(typ(p6)==10)&&(lgef(p6)>3);
  1011.       if(f) goto tchi;
  1012.       p1=factor(p5);p2=(GEN)p1[1];l=lg(p2)-1;
  1013.       switch(l)
  1014.         {
  1015.         case 1: f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
  1016.           if(f) {y[2]=un;y[1]=lstoi(12);return y;}
  1017.           else {y[2]=lneg(gun);y[1]=lstoi(24);return y;}
  1018.         case 2: avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(8);return y;
  1019.         case 3: l2=lgef(p2[1])-3;
  1020.           if(l2==2) {avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(4);return y;}
  1021.           else {avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(4);return y;}
  1022.         default: err(talker,"bug1 in galois");
  1023.         }
  1024.     case 5:
  1025.       do
  1026.         {
  1027.           do
  1028.         {
  1029.           p1=roots(x,prec);z=cgetg(7,17);
  1030.           for(l=1;l<=5;l++)
  1031.             {
  1032.               p2=(l==1)?p1:transroot(p1,1,l);
  1033.               p3=gzero;k=0;for(i=1;i<=5;i++)
  1034.             {
  1035.               p5=gadd(gmul(p2[ind5[k]],p2[ind5[k+1]]),gmul(p2[ind5[k+2]],p2[ind5[k+3]]));
  1036.               p3=gadd(p3,gmul(gsqr(p2[i]),p5));k+=4;
  1037.             }
  1038.               z[l]=lrndtoi(greal(p3),&e);
  1039.               p4=(l==1)?gsub(polx[0],p3):gmul(p4,gsub(polx[0],p3));
  1040.             }
  1041.           p2=transroot(p1,2,5);
  1042.           p3=gzero;k=0;for(i=1;i<=5;i++)
  1043.             {
  1044.               p5=gadd(gmul(p2[ind5[k]],p2[ind5[k+1]]),gmul(p2[ind5[k+2]],p2[ind5[k+3]]));
  1045.               p3=gadd(p3,gmul(gsqr(p2[i]),p5));k+=4;
  1046.             }
  1047.           z[6]=lrndtoi(greal(p3),&e);
  1048.           p4=gmul(p4,gsub(polx[0],p3));
  1049.           p5=grndtoi(greal(p4),&e);
  1050.           e=max(e,gexpo(gimag(p4)));
  1051.           if(e> -10) prec=(prec<<2)-2;
  1052.         }
  1053.           while(e> -10);
  1054.           p6=ggcd(p5,deriv(p5,0));
  1055.           f=(typ(p6)==10)&&(lgef(p6)>3);
  1056.           if(f) goto tchi;
  1057.           p3=factor(p5);l=lg(p3[1])-1;
  1058.           f=carreparfait(discsr(x));
  1059.           if(l==1)
  1060.         {
  1061.           avma=av;y=cgetg(4,17);y[3]=un;
  1062.           if(f) {y[2]=un;y[1]=lstoi(60);return y;}
  1063.           else {y[2]=lneg(gun);y[1]=lstoi(120);return y;}
  1064.         }
  1065.           else
  1066.         {
  1067.           if(f) 
  1068.             {
  1069.               l=1;while((l<=6)&&(!gcmp0(poleval(p5,z[l])))) l++;
  1070.               if(l>6) err(talker,"bug4 in galois");
  1071.               p2=(l==6)?transroot(p1,2,5):transroot(p1,1,l);
  1072.               p3=gzero;
  1073.               for(i=1;i<=5;i++)
  1074.             {
  1075.               j=(i%5)+1;
  1076.               p3=gadd(p3,gmul(gmul(p2[i],p2[j]),gsub(p2[j],p2[i])));
  1077.             }
  1078.               p5=gmul(p3,p3);p4=grndtoi(greal(p5),&e1);
  1079.               e1=max(e1,gexpo(gimag(p5)));
  1080.               if(e1<= -10)
  1081.             {
  1082.               if(gcmp0(p4)) goto tchi;
  1083.               f=carreparfait(p4);
  1084.               avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(f?5:10);return y;
  1085.             }
  1086.               else prec=(prec<<2)-2;
  1087.             }
  1088.           else
  1089.             {
  1090.               avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(20);return y;
  1091.             }
  1092.         }
  1093.         }
  1094.       while(e1> -10);
  1095.     case 6: 
  1096.       do
  1097.         {
  1098.           do
  1099.         {
  1100.           p1=roots(x,prec);
  1101.           for(l=1;l<=6;l++)
  1102.             {
  1103.               p2=(l==1)?p1:transroot(p1,1,l);
  1104.               p3=gzero;k=0;for(i=1;i<=5;i++) for(j=i+1;j<=6;j++)
  1105.             {
  1106.               p5=gadd(gmul(p2[ind6[k]],p2[ind6[k+1]]),gmul(p2[ind6[k+2]],p2[ind6[k+3]]));
  1107.               p3=gadd(p3,gmul(gsqr(gmul(p2[i],p2[j])),p5));k+=4;
  1108.             }
  1109.               p4=(l==1)?gsub(polx[0],p3):gmul(p4,gsub(polx[0],p3));
  1110.             }
  1111.           p5=grndtoi(greal(p4),&e);
  1112.           e=max(e,gexpo(gimag(p4)));
  1113.           if(e> -10) prec=(prec<<2)-2;
  1114.         }
  1115.           while(e> -10);
  1116.           p6=ggcd(p5,deriv(p5,0));
  1117.           f=(typ(p6)==10)&&(lgef(p6)>3);
  1118.           if(f) goto tchi;
  1119.           p3=factor(p5);p2=(GEN)p3[1];l=lg(p2)-1;
  1120.           switch(l)
  1121.         {
  1122.         case 1:    p3=gadd(gmul(gmul(p1[1],p1[2]),p1[3]),gmul(gmul(p1[4],p1[5]),p1[6]));
  1123.           p4=gsub(polx[0],p3);
  1124.           for(i=1;i<=3;i++)
  1125.             for(j=4;j<=6;j++)
  1126.               {
  1127.             p2=transroot(p1,i,j);
  1128.             p3=gadd(gmul(gmul(p2[1],p2[2]),p2[3]),gmul(gmul(p2[4],p2[5]),p2[6]));
  1129.             p4=gmul(p4,gsub(polx[0],p3));
  1130.               }
  1131.           p5=grndtoi(greal(p4),&e1);
  1132.           e1=max(e1,gexpo(gimag(p4)));
  1133.           if(e1<= 10)
  1134.             {
  1135.               p6=ggcd(p5,deriv(p5,0));
  1136.               f=(typ(p6)==10)&&(lgef(p6)>3);
  1137.               if(f) goto tchi;
  1138.               p3=factor(p5);p2=(GEN)p3[1];l=lg(p2)-1;f=carreparfait(discsr(x));
  1139.               avma=av;y=cgetg(4,17);y[3]=un;
  1140.               if(l==1)
  1141.             {
  1142.               if(f) {y[2]=un;y[1]=lstoi(360);return y;}
  1143.               else {y[2]=lneg(un);y[1]=lstoi(720);return y;}
  1144.             }
  1145.               else
  1146.             {
  1147.               if(f) {y[2]=un;y[1]=lstoi(36);return y;}
  1148.               else {y[2]=lneg(un);y[1]=lstoi(72);return y;}
  1149.             }
  1150.             }
  1151.           else prec=(prec<<2)-2;
  1152.           break;
  1153.           
  1154.         case 2: l2=lgef(p2[1])-3;if(l2>3) l2=6-l2;
  1155.           switch(l2)
  1156.             {
  1157.             case 1: f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
  1158.               if(f) {y[2]=un;y[1]=lstoi(60);return y;}
  1159.               else {y[2]=lneg(un);y[1]=lstoi(120);return y;}
  1160.             case 2: f=carreparfait(discsr(x));
  1161.               if(f) {avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(24);return y;}
  1162.               else
  1163.             {
  1164.               p3=(lgef(p2[1])==5)?(GEN)p2[2]:(GEN)p2[1];
  1165.               f=carreparfait(discsr(p3));avma=av;y=cgetg(4,17);y[2]=lneg(gun);
  1166.               if(f) {y[1]=lstoi(24);y[3]=deux;return y;}
  1167.               else {y[1]=lstoi(48);y[3]=un;return y;}
  1168.             }
  1169.             case 3: f=carreparfait(discsr(p2[1]))||carreparfait(discsr(p2[2]));
  1170.               avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(f?18:36);
  1171.               return y;
  1172.             }
  1173.         case 3: for(l2=1;l2<=3;l2++) if(lgef(p2[l2])>=6) p3=(GEN)p2[l2];
  1174.           if(lgef(p3)==6)
  1175.             {
  1176.               f=carreparfait(discsr(p3));avma=av;y=cgetg(4,17);y[2]=lneg(gun);
  1177.               y[3]=un;y[1]=f?lstoi(6):lstoi(12);return y;
  1178.             }
  1179.           else
  1180.             {
  1181.               f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
  1182.               if(f) {y[2]=un;y[1]=lstoi(12);return y;}
  1183.               else {y[2]=lneg(un);y[1]=lstoi(24);return y;}
  1184.             }
  1185.         case 4: avma=av;y=cgetg(4,17);y[1]=lstoi(6);y[2]=lneg(gun);
  1186.           y[3]=deux;return y;
  1187.         default: err(talker,"bug3 in galois");
  1188.         }
  1189.         }
  1190.       while(e1> -10);
  1191.       
  1192.     case 7: 
  1193.       do
  1194.         {
  1195.           p1=roots(x,prec);p4=gun;
  1196.           for(i=1;i<=5;i++)
  1197.         for(j=i+1;j<=6;j++)
  1198.           for(k=j+1;k<=7;k++)
  1199.             p4=gmul(p4,gsub(polx[0],gadd(gadd(p1[i],p1[j]),p1[k])));
  1200.           p5=grndtoi(greal(p4),&e);e=max(e,gexpo(gimag(p4)));
  1201.           if(e> -10) prec=(prec<<2)-2;
  1202.         }
  1203.       while(e> -10);
  1204.       p6=ggcd(p5,deriv(p5,0));
  1205.       f=(typ(p6)==10)&&(lgef(p6)>3);
  1206.       if(f) goto tchi;
  1207.       p1=factor(p5);p2=(GEN)p1[1];l=lg(p2)-1;
  1208.       switch(l)
  1209.         {
  1210.         case 1: f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
  1211.           if(f) {y[2]=un;y[1]=lstoi(2520);return y;}
  1212.           else {y[2]=lneg(gun);y[1]=lstoi(5040);return y;}
  1213.         case 2: f=lgef(p2[1])-3;avma=av;y=cgetg(4,17);y[3]=un;
  1214.           if((f==7)||(f==28)) {y[2]=un;y[1]=lstoi(168);return y;}
  1215.           else {y[2]=lneg(gun);y[1]=lstoi(42);return y;}
  1216.         case 3: avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(21);return y;
  1217.         case 4: avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(14);return y;
  1218.         case 5: avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(7);return y;
  1219.         default: err(talker,"bug2 in galois");
  1220.         }
  1221.     }
  1222.     tchi:
  1223.       avma=av1;x=tchirnhausen(x1);
  1224.     }
  1225. }
  1226.  
  1227. GEN initalg(x,prec)
  1228.      GEN x;
  1229.      long prec;
  1230. {
  1231.   GEN y,p1,p2,p3,p4,p5,p6,p7,fieldd,polr,ptrace,dx,adx,polmax,s,a,dxn,adxn,sn,phimax,ind;
  1232.   long tx=typ(x),n=lgef(x)-3,i,j,av=avma,av1,v,k,r1,imax,numb,flc,tetpil;
  1233.  
  1234.   if((tx!=10)||(n<=0)) err(poltyper);
  1235.   p1=content(x);p1=gcmp1(p1) ? x: gdiv(x,content(x));
  1236.   for(k=2;k<=n+2;k++) if(typ(p1[k])!=1) err(impl,"general algebraic extension");
  1237.   if(gcmp1(p2=(GEN)p1[n+2])) x=p1;
  1238.   else
  1239.     {
  1240.       x=cgetg(n+3,10);x[1]=p1[1];x[n+2]=un;x[n+1]=p1[n+1];p3=p2;
  1241.       for(k=n;k>=2;k--)
  1242.     {
  1243.       x[k]=lmulii(p3,p1[k]);
  1244.       if(k>2) p3=mulii(p2,p3);
  1245.     }
  1246.     }
  1247.   p1=factor(x);
  1248.   if(lgef(coeff(p1,1,1))!=n+3) err(redper1);
  1249.   r1=sturm(x);p4=allbase(x,0,&fieldd);
  1250.   if(r1<n)
  1251.     {
  1252.       polr=roots(x,prec);p3=cgetg(n+1,19);
  1253.       for(i=1;i<=n;i++)
  1254.     {
  1255.       p1=cgetg(n+1,18);p3[i]=(long)p1;
  1256.       for(j=1;j<=n;j++)
  1257.         p1[j]=lsubst(p4[i],varn(p4[i]),polr[j]);
  1258.     }
  1259.       p2=greal(gmul(gconj(gtrans(p3)),p3));
  1260.     }
  1261.   else
  1262.     {
  1263.       ptrace=cgetg(n+1,17);ptrace[1]=lstoi(n);
  1264.       for(k=1;k<n;k++) 
  1265.     {
  1266.       p3=gmulsg(k,x[n-k+2]);
  1267.       for(i=1;i<k;i++) p3=gadd(p3,gmul(x[n-i+2],ptrace[k-i+1]));
  1268.       ptrace[k+1]=lneg(p3);
  1269.     }
  1270.       p2=cgetg(n+1,19);
  1271.       for(i=1;i<=n;i++)
  1272.     {
  1273.       p1=cgetg(n+1,18);p2[i]=(long)p1;
  1274.       for(j=1;j<i;j++) p1[j]=lcopy(coeff(p2,i,j));
  1275.       for(j=i;j<=n;j++)
  1276.         {
  1277.           p5=gres(gmul(p4[i],p4[j]),x);p6=gzero;
  1278.           for(k=0;k<=lgef(p5)-3;k++) p6=gadd(p6,gmul(p5[k+2],ptrace[k+1]));
  1279.           p1[j]=(long)p6;
  1280.         }
  1281.     }
  1282.     }
  1283.   p1=lllgram(p2,prec);v=varn(x);dx=discsr(x);adx=absi(dx);polmax=x;imax=0;
  1284.   if(r1<n) for(s=gzero,i=1;i<=n;i++) s=gadd(s,gnorm(polr[i]));
  1285.   else s=gsub(gsqr(x[n+1]),gmul2n(x[n],1));
  1286.   a=cgetg(n+1,18);for(i=1;i<=n;i++) a[i]=lmul(p4,p1[i]);
  1287.   for(numb=0,i=1;i<=n;i++)
  1288.     {
  1289.       av1=avma;p3=gmodulcp(a[i],x);p7=content(p3[2]);
  1290.       if(gcmp1(p7)) p3=caract(p3,v);
  1291.       else
  1292.     {
  1293.       p3=caract(gdiv(p3,p7),v);
  1294.       p3=gmul(gpuigs(p7,lgef(p3)-3),gsubst(p3,v,gdiv(polx[v],p7)));
  1295.     }
  1296.       p5=ggcd(deriv(p3,v),p3);
  1297.       if(lgef(p5)==3)
  1298.     {
  1299.       dxn=discsr(p3);adxn=absi(dxn);flc=gcmp(adxn,adx);numb++;
  1300.       if(flc<=0)
  1301.         {
  1302.           if(r1<n) for(sn=gzero,j=1;j<=n;j++) sn=gadd(sn,gnorm(poleval(a[i],polr[j])));
  1303.           else sn=gsub(gsqr(p3[n+1]),gmul2n(p3[n],1));
  1304.           if(flc<0) {dx=dxn;adx=adxn;s=sn;polmax=p3;imax=i;}
  1305.           else 
  1306.         {
  1307.           flc=gcmp(sn,s);
  1308.           if((flc<0)||((!flc)&&(gpolcomp(p3,polmax)<0)))
  1309.             {dx=dxn;adx=adxn;s=sn;polmax=p3;imax=i;}
  1310.         }
  1311.         }
  1312.     }
  1313.     }
  1314.   if(!numb) err(talker,"you have found a counterexample to a conjecture,\nplease send us the polynomial as soon as possible");
  1315.   phimax=imax?(GEN)a[imax]:polx[v];
  1316.   j=n+1;
  1317.   while((j>=2)&&(!signe(polmax[j]))) j-=2;
  1318.   if((j>=2)&&(signe(polmax[j])>0))
  1319.     {
  1320.       for(;j>=2;j-=2) setsigne(polmax[j],-signe(polmax[j]));
  1321.       phimax=gneg(phimax);
  1322.     }
  1323.   p2=lift(gsubst(p4,v,polymodrecip(gmodulcp(phimax,x))));
  1324.   p3=cgetg(n+1,19);
  1325.   for(j=1;j<=n;j++)
  1326.     {
  1327.       p4=cgetg(n+1,18);p3[j]=(long)p4;
  1328.       for(i=1;i<=n;i++) p4[i]=(long)truecoeff(p2[j],i-1);
  1329.     }
  1330.   p4=denom(p3);
  1331.   p2=gdiv(hnf(gmul(p4,p3)),p4);p3=cgetg(n+1,17);
  1332.   for(j=1;j<=n;j++) 
  1333.     {
  1334.       p1=gzero;for(i=n;i;i--) p1=gadd(coeff(p2,i,j),gmul(p1,polx[v]));
  1335.       p3[j]=(long)p1;
  1336.     }
  1337.   if(!carrecomplet(divii(dx,fieldd),&ind)) err(talker,"bug1 in initalg");
  1338.   tetpil=avma;y=cgetg(8,17);y[1]=lcopy(polmax);p1=cgetg(3,17);
  1339.   p1[1]=lstoi(r1);p1[2]=lstoi((n-r1)>>1);y[2]=(long)p1;
  1340.   y[3]=lcopy(fieldd);y[4]=lcopy(ind);y[5]=lcopy(s);
  1341.   if(r1<n) 
  1342.     {
  1343.       p1=cgetg(n+1,18);
  1344.       for(i=1;i<=n;i++) p1[i]=(long)poleval(phimax,polr[i]);
  1345.       y[6]=(long)p1;
  1346.     }
  1347.   else y[6]=lreal(roots(polmax,prec));
  1348.   y[7]=lcopy(p3);
  1349.   return gerepile(av,tetpil,y);
  1350. }
  1351.  
  1352. int gpolcomp(p1,p2)
  1353.      GEN p1,p2;
  1354. {
  1355.   int d,j;
  1356.  
  1357.   d=lgef(p1)-3;if((lgef(p2)-3)!=d) err(talker,"bug: different degrees in gpolcomp");
  1358.   j=d+1;while((j>=2)&&gegal(absi(p1[j]),absi(p2[j]))) j--;
  1359.   if(j==1) return 0;
  1360.   return gcmp(absi(p1[j]),absi(p2[j]));
  1361. }
  1362.  
  1363. GEN galoisconj2(x,prec)
  1364.      GEN x;
  1365.      long prec;
  1366. {
  1367.   long av=avma,tetpil,i,j,n,v;
  1368.   GEN y,w,polr,p1,p2,b,di;
  1369.  
  1370.   if(typ(x)!=10) err(galer1);
  1371.   n=lgef(x)-3;if(n<=0) err(galer1);
  1372.   v=varn(x);p1=factor(x);
  1373.   if((lg(p1[1])>2)||(!gcmp1(coeff(p1,1,2)))) err(galcer1);
  1374.   polr=roots(x,prec);p1=(GEN)polr[1];b=allbase(x,0,&di);
  1375.   w=cgetg(n+1,17);w[1]=un;for(i=2;i<=n;i++) w[i]=(long)gsubst(b[i],v,p1);
  1376.   y=cgetg(n+1,17);y[1]=(long)polx[v];
  1377.   for(i=2;i<=n;i++)
  1378.     {
  1379.       p1=lindep(concat(w,polr[i]),prec);
  1380.       if(gcmp0(p1[n+1])) y[i]=zero;
  1381.       else
  1382.     {
  1383.       p2=gzero;
  1384.       for(j=1;j<=n;j++) p2=gadd(p2,gmul(p1[j],b[j]));
  1385.       p2=gdiv(p2,gneg(p1[n+1]));
  1386.       if(gcmp0(gmod(gsubst(x,v,p2),x))) y[i]=(long)p2;else y[i]=zero;
  1387.     }
  1388.     }
  1389.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1390. }
  1391.  
  1392. GEN galoisconj(x,prec)
  1393.      GEN x;
  1394.      long prec;
  1395. {
  1396.   long av=avma,tetpil,i,j,n,v;
  1397.   GEN y,w,polr,p1,p2;
  1398.  
  1399.   if(typ(x)!=10) err(galer1);
  1400.   n=lgef(x)-3;if(n<=0) err(galer1);
  1401.   v=varn(x);p1=factor(x);
  1402.   if((lg(p1[1])>2)||(!gcmp1(coeff(p1,1,2)))) err(galcer1);
  1403.   polr=roots(x,prec);p1=(GEN)polr[1];
  1404.   w=cgetg(n+1,17);w[1]=un;for(i=2;i<=n;i++) w[i]=lmul(p1,w[i-1]);
  1405.   y=cgetg(n+1,17);y[1]=(long)polx[v];
  1406.   for(i=2;i<=n;i++)
  1407.     {
  1408.       p1=lindep(concat(w,polr[i]),prec);
  1409.       if(gcmp0(p1[n+1])) y[i]=zero;
  1410.       else
  1411.     {
  1412.       p2=gzero;
  1413.       for(j=n;j;j--) p2=gadd(p1[j],gmul(p2,polx[v]));
  1414.       p2=gdiv(p2,gneg(p1[n+1]));
  1415.       if(gcmp0(gmod(gsubst(x,v,p2),x))) y[i]=(long)p2;else y[i]=zero;
  1416.     }
  1417.     }
  1418.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1419. }
  1420.  
  1421.  
  1422.  
  1423.  
  1424.  
  1425.  
  1426.  
  1427.  
  1428.  
  1429.  
  1430.